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Abstract 


A new parabolized Navier-Stokes (PNS) code has been 
developed to compute the hypersonic, viscous, chemically 
reacting flow fields around three-dimensional bodies. The' 
flow medium is assumed to be a multicomponent mixture of 
thermally perfect but calorically imperfect gases. The new 
PNS code solves the gas dynamic and species conservation 
equations in a coupled manner using a noniterative, im- 
plicit, approximately-factored, finite-difference algorithm. 
The space-marching method is made well-posed by special 
treatment of the streamwise pressure gradient term. The 
code has been used to compute hyperaonic laminar flow of 
chemically reacting air over cones at angles of attack. The 
results of the computations are compared with the results of 
reacting boundary-layer computations and show excellent 
agreement. 
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NomcnsJsturs 

reference area, (m 1 ) 

frozen speed of sound 

total drag coefficient 

skin-friction coefficient 

pressure coefficient 

frozen specific heat of the mixture 

mass fraction of species a 

kinematic binary diffusion coefficient 

static enthalpy 

total enthalpy of the mixture 

enthalpy of formation of species * 

backward reaction rate constant for 

the mth reaction 

forward reaction rate constant for 
the mth reaction 
reference length, (m) 
binary Lewis number 
diffusion mass flux vector 
number of reactions 
frozen Mach number 
molecular mass 
number of species 
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Avogadro number, 6.022169 x 10 1 '’ ( kmol' *) 
electron number density, (m~ 3 ) 
number of reactants (including catalytic 
third bodies) 

static pressure of the mixture 
heat flux vector 

universal gas constant, 8314.34 Jftkmal K] 
Reynolds number based on L‘ 

Stanton number 
time 

static temperature of the mixture 
Cartesian components of the mass-averaged 
velocity 

magnitude of the mass-averaged velocity 

mass production rate of species s 

mole fraction of species s 

Cartesian coordinates (physical space) 

third body efficiency (relative to argon) 

of the rth catalytic body 

angle of attack 

mole-mass ratio of species a 

cone half angle 

thermal conductivity 

molecular viscoaity 

generalised coordinates 

mass density 

safety factor 

shear stress 

finite-difference indices 
indices denoting species 
wail 

partial derivative w.r.t. z,y,z 
partial derivative w.r.t. £,»j.c 
freestream 


Superscripts 
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e : 

i : 
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*,y.* 
£.*?. f 


dimensional quantity 
chemical quantity 
inviscid quantity 
transpose 
viscous quantity 
Cartesian components 
transformed components 


Intro duction 

Accurate numerical simulations of the acrothcrmody- 
namic environments around the recently proposed hyper- 
sonic vehicles’ * must account for the various complex re- 
laxation phenomena such as vibrational excitation, chem- 
ical reactions, ionization, tie., which occur in the shock 
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layers. The numerical simulation of complete thermochem- 
leal nonequilibrium has been successfully accomplished by 
Park*' 4 for one-dimensional flows. The extension of this 
analysis to higher dimensional flows is quite complicated 
and requires further research. However, as a first step, the 
problem can be simplified by assuming the flow field to 
be in thermal equilibrium but not in chemical equilibrium. 
The numerical simulation of viscous, chemically reacting, 
external flows around three-dimensional configurations is 
the focus of attention of several investigations today. 

The numerical methods currently employed fall into two 
main categories, (i) time-marching methods and (ii) space- 
marching methods. In the time-marching methods,*"* time- 
asymptotic solutions to the Navier-Stokes equations are 
; computed. These methods are accurate but require a sub- 
! stantial amount of computer time. Space-marching meth- 
ods, on the other hand, require much less computer time 
and provide accurate solutions in cases where they are ap- 
plicable. In the latter category, the viscous shock-layer 
(VSL) equations 8,10 have been widely used to compute 
three-dimensional, viscous, chemically reacting flows. These 
equations are uniformly valid in the shock layer but fail 
in the presence of crossflow separation. This deficiency is 
overcome through the use of the parabolized Navier-Stokes 
(PNS) equations. Bhutta and Lewis 11 were one of the first 
to use the PNS equations to compute chemical nonequilib- 
rium flow fields. They solved the chemistry and gas dynam- 
ics separately and used an iterative approach to couple the 
two. In contrast to their approach, Prabhu ct al. li devel- 
oped a chemical nonequilibrium PNS code in which the gaa 
dynamics and chemistry were solved simultaneously in a 
noniterative manner. This two-dimensional /axisymmetrk 
code was successfully applied to reacting flows over wedgee 
and cones and was found to be very competitive in terms 
of both efficiency as well as computing time. Having es- 
tablished the feasibility of computing reacting ftowa with a 
coupled approach, a new PNS code for three-dimenalonal 
bodies has been developed in the present study. 

The present paper describes in detail the development 
of this new PNS code which can compute chemically react- 
ing How fields around three-dimensional bodies. As in the 
previous study 'L the gas dynamic and species conservation 
equations are solved in a coupled manner. A noniterative, 
implicit, approximately-factored, finite-difference method 
is used to solve the PNS equations. The flow medium con- 
sidered in the present computations is air contisting of aix 
species and electrons. The electrons are eliminated from 
the species set using the principle of conservation of charge. 
The source terms are treated in a partially implicit man- 
ner in the present formulation. The code has been used to 
compute a number of chemically reacting flow fields around 
cones at angles of attack. The results of the computa- 
tions are compared with those of a reacting boundary-Uyer 
code. 13 

Govern! ng E quations 

The equations governing the steady, three-dimensional 
laminar flow of a multicomponent gaa have been obtained 
from the equations in Ref. 14 by (i) neglecting the time- 
derivative terms, (ii) assuming the flow to be in thermal 
equilibrium, (iii) neglecting radiation, and by (iv) assum- 
ing mass diffusion to be binary and due to concentration 


gradients only. These equations can be written in nondi- 
mensional, strong conservation-law form in Cartesian coor- 
dinates for an n-component system as: 

E’, + Fi-rGi= (1) 

The (n + 4)-component vector of conservation variables 
Q is chosen as 

Q = {p,pu,pt),pw,p/f,p<i,pca,...,pc R _i} r (2) 

E\ F\ and G* are the (n + 4)-<omponent inviscid flux 
vectors. E", F“,and G“ are (n - 4)-component viscous flux 
vectors. W e ir the vector of chemical source terms. All 
these vector* are functions of the elements of the vector Q 
and their spatial derivatives. Explicit expressions for ele- 
ments of the flux and source vectors are given in Appendix 
A. Note that the vector of dependent variables contains 
both the fluid dynamic and chemistry variables and that 
only (n - 1) of the n species continuity equations are re- 
quired because the mass fractions sum to unity. Thus, the 
nth species continuity equation is replaced by the follow ing 
algebraic equation 
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e„ = i y 


( 3 ) 
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where Eq. 4 is the equation of state for perfect gases, Eq. 
S ia the definition of the mixture total enthalpy and Eq. 6 
is the definition of the mixture molecular mass. 

The following nondimensionalisation has been employed 
in the present formulation 
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The other nondimensional quantities appearing in the 
equations are 
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Iii order to close the system of governing equations, the 
liirriiiodynamir and transport properties of the constituent 
gases and the mixture are required. This is discussed in the 
next subsection. 

Gas Model. Thermodynamic and Transport Pr operties 
(a) Gas Model 

The chemical model used in the present calculations is 
air consisting of molecular oxygen (Oj), atomic oxygen (0), 
molecular nitrogen (A'j), nitric oxide (NO), nitric oxide ion 
(NO~). atomic nitrogen (A'), and electrons («*). These 
species are indexed s 1 - 6 in the order shown with the 
electrons being treated as a special case. The following 
reactions are considered between the constituent species 


(1) 

O j + Mi 

s* 20 

+ Mi 

(2) 

Nj + Mt 

- 2N 

■+ M 2 

(3) 

N 2 + N 

«=* 2N 

+ N 

(4) 

NO + Ms 

** N 

+ 0 + Ms 

(3) 

NO + 0 

~ o 2 

+ N 

(6) 

n 3 +o 

~ NO 

+ N 

(7) 

N +0 

« NO* 

+ e~ 


where M 1( AS], and Ms are catalytic third bodies. The 
above model has six species (n — 6), seven reactions (m = 
7) and ten reactants (n, = 10) including electrons. The 
reaction rates have been obtained from Blottner. 18 

(b) Specific Heat and En t halpy 

The enthalpies ( J/kg ) and specific heats (J/kg-K) of 
the species are obtained from the following relations 

*; = rc lt ,(r) + c 

c;.. = c 2 ,,(n 

Tables of Ci t# and C 2i , as functions of 7" (K) are obtained 
from Ref. 13. Cubic spline interpolation is used in these 
tables. The enthalpy and frozen specific heat of the mixture 
are given by the following expressions 


(io) 

*= i 


fj 


d/r 1 

dT- j. 


v- dh: 

r L r ' 


»- 1 


dT 




r- 1 


where the subscripts on the differenlistion indicate that the 
composition of the mixture is frozen locally. 


(c) Viscosity and Thermal Conductivity 

The viscosity (A'-s/m 3 ) of species a is calculated from 
curve fils developed in Ref. 13. These curve fits are of the 
form 


t‘] 0. 1 exp ; ( .4 , log, T - /f,)log, T ♦ (’,] (12) 


where A«, B„ and C, are constants. The thermal conduc- 
tivity ( W l m-K ) of species a is computed using Eucken’s 
semiempiriral formula 
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(13) 


The viscosity and thermal conductivity of the mixture 
are calculated using Wilke's semiempirica! mixing rule 15 
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(16) 

Wilke's mixing rule is considered adequate for weakly ion- 
ising flows. 

( d) Diffusion Coefficient 

The binary Lewis numbers for all the species are as- 
sumed to be the same constant Ce. The kinematic binary 
diffusion coefficient D~ (m 3, .«) is then computed from the 
definition 

k‘ Lt 


D ■ = 




(17) 


The dimensional thermodynamic and transport properties 
are nondimensionalised using Eq. 7 prior to their use in 
the code. 


Consider a generalized transformation of the spatial co- 
ordinates as shewn below 


£=£(*,¥.*). q = rj(x,y.z), f = f(*,y,z) ( 18 ) 


where £ is the “streamwise" coordinate. r\ is the ‘’normal’' 
coordinate and f is the “circumferential" coordinate. Using 
this transformation, the governing equations ran bp recast 
into the strong conservation form given by 

+ f; i g; * - (ev - r„ - g?) - w (19) 

where E', f 6', £\ l >v , G w , and W c , are the transformed 
inviacid and viscous fluxes and chemical source terms, re- 
spectively. Their forms are given in Appendix A. 

The equation set, Eq. 19, is simplified by making the 
thin-layer approximation, ».e., the viscous snd diffusion ef- 
fects in the streamwise and meridional directions are as- 
sumed to be negligibly small compared to those in the nor- 
mal direction. Therefore, the set of equations is simplified 
to 

+ £*„ iG!= 4 W c (20) 

tit 

where t contains derivatives with respect to the q coor- 
dinate only. 

Streamwise Pressure Gradient 

In its present form, Eq. 20, is hyperbolic-elliptic in the 
streamwise or marching direction £ and consequently the 
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space-marching method of solution is ill-posed. The tech- 
nique used to overcome this problem tvs a originally pro- 
posed by Vigneron et ol. 10 for ideal gases and later extended 
to chemically reacting flows by Prabhu et a/. 13 Using this 
technique it was found that the governing equations are 
hyperbolic-parabolic if and only if (i) there is no axial flow 
separation in the solution domain, (ii) the local frocen Mach 
number is greater than unity in the inviscid part of the flow 
held, and (iii) only a fraction u? of the streamwise pressure 
gradient dp/di is retained in the subsonic part of the flow 
field. The magnitude of w depends on the local frosen Mach 
number and is determined from the expression 


u = min 


jl.oM/ 5 1 -4 v(A if - l) | (21) 


where o (0.8 < o < 0.9) is a factor of safety and 

0. 


X - 


Mj 


MC„, 
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(23) 

(24) 


' (1 A) 

Introducing w into Eq. 20. the streamwise inviscid flux 
can be split as follows 


E’ - E' 4 £'' 


(26) 


where 


£•'" = (1 - w)p{0, ^,0,0,0 o} r (26) 


Using Eq. 26, Eq. 20 can be rewritten as 


E< + K + 6} = - *>«" (27) 

Further, the last term of Eq. 27 is normally neglected 
in the subsonic region. Hence the final set of equations that 
is solved is 

£•' + *• + G; = + W' ( 28 ) 


Numerical Solution of PNS Equat ions 
Finit e-difference Algorithm 

The numerical algorithm used to solve the system of 
equations, Eq. 28, is an adaptation of the one developed 
by Tannehill et of.’ 7 The factored algorithm is implemented 
as the following sequence of steps 
(a) Normal Sweep: k =■ 2^3 



(b) Circumferential Sweep: j = 1,2, ■ ■ ■ , N J 



- A(Af >tj + Af— (B iJtJ - 


- (A,,*j - AfA^jAQ,.*., (30) 

(c) Update: 


Q = Qijkj + AQ iJkj - (31) 


where the Jacobian matrices are 


A = 
M = 


dfe’ 1 
dQ • 

&y 

dQ * 


B = 


A e = 


dt ( 

dQ' 

dW e 

aq 


ae* 

dQ 


(32) 


The subscripts i,k,j are indices associated with the di- 
rections (,q,f, respectively. NJ and KK are the total 
number of grid points in the c and q directions, respec- 
tively. A( is the marching stepsise. The derivatives d/dq 
and d/df are replaced by conventional three-point central- 
difference operators. The algorithm is first-order accurate 
in the { direction and second-order accurate in the q and f 
directions. FVeestream fluxes are subtracted from the invis- 
cid fluxes in order to preserve freest re am. Fourth-difference 
explicit and second-difference implicit smoothing operators 
are also added to the factored operators. These have not 
been shown above but their forms can be found in Ref. 12. 

The Jacobian matrices A, B, C, M, and A e represent 
the linearisation of the fluxes and source terms and the 
elements of these matrices are given in Appendices B, C, 
and D. In the linearisation of the viscous flux, the transport 
properties are assumed to be locally constant. The caret 
above the symbols for the fluxes, Jacobians, and source 
terms signifies that the geometry is not linearized along 
with the flow variables. 

The left hand sides of Eqs. 29 and 30 correspond to 
a block-tridiagonal system of equations. The blocks are 
square matrices of order (n + 4). For the six-species air 
model considered in the present calculations, the blocks are 
square matrices of order 10. The block-tridiagonal solver 
for the 10 x 10 blocks is developed along the same lines as 
the one by Steger* 8 for 5 x 5 blocks. 


Boundary Conditions 

In the present work, only flows without yaw are consid- 
ered. Therefore, at every streamwise station the computa- 
tional domain is bounded by (i) the outer boundary which 
is taken to be the freest ream, (ii) the inner boundary which 
is taken to be the wall, and (iii) the pitch plane of symme- 
try. Aqy discontinuities in. the fiowfield are “captured" as 
a part of the solution. 

At the pitch plane of symmetry reflection boundary con- 
ditions are imposed and thus, flow symmetry is maintained. 

The following nondimenaional boundary conditions are 
imposed implicitly at the wall 
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(a) u ~ 0* v - 0, w - 0, fa (p) - 0 

(b) T - T w (isothermal wall) or fa (7') - 0 (adiabatic wail) 

(c) t, = c,„ (catalytic wall) or fa (e.) = 0 (noncaulytic 
wall), a = 1,2,..., n- 1 

The nondimensional boundary conditions at the outer 
boundary are 

(a) u = coso, u = 0, u> = sin a 

(b) T = 1, /> = 1, e, = e tm , a = 1,2 n - 1 

Initial C ond itions 

The PNS equations require initial conditions in addition 
to the boundary conditions. The usual procedure is to use 
an initial data surface generated by a full Navier-Stokee 
code. For conical or pointed bodies, however, the code 
generates it own starting solution. This starting solution is 
generated iteratively using a “step back" procedure. ,2,l# 

Decoding 

The primitive variables p,u,v,w,H,e i,cj,...,c,_| at 
station t + 1 are easily obtained from the elements of Q. 
The mass fraction of the nth species is computed using 
Cq. 3 and the static enthalpy of the mixture is computed 
using Gq. 5. For a given species distribution t j, c 3 , ..., 
c„ and mixture enthalpy h, the temperature is iteratively 
determined using the following algorithm 12 


rp'k - ♦ I _ 


ZlVc'CrAT'*) 


(S3) 


where k is the index of iteration. The iterations are contin- 
ued until 

|r* +1 -r* !«« ( 34 ) 

where t is a small positive quantity. Once the temperature 
is determined, the thermodynamic and transport proper- 
ties are easily computed using the expressions given in the 
previous sections. 


Grid Generation 


An algebraic grid generation procedure is used in the 
present calculations. In this procedure, the point on the 
body surface and the point on the outer boundary are con- 
nected by a straight line and the grid points are distributed 
on this line using the following stretching function 


(*-0 
(SK - 1) 


j'-” -(0-1)'-') 

(33) 

+ (0 - 1)'-' j 

k -- 1,2 ,...,NK 

(36) 


where 0 (0 > 1) is the stretching parameter and NK is the 
total number of points on the grid line. Note that *(0) = 0 
and s(l) = ] and the points are clustered close to the wall 
for values of 0 close to 1. Such clustering is necessary for 
good resolution of the subsonic viscous layer. 

The coordinates of the grid points are then obtained 
from 
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* = x w + r»is(ij)4({) 

y - yu. + ( 3 ~) 

* = tw + nsa{ti)6(Q 

where «i, n 3 , and r»j are the direction cosines of the unit 
vector along the grid tine and 6 is the linear distance of the 
outer boundary from the body surface. 

The metrics are then computed using first-order accu- 
rate one-sided differences in the £ direction and second- 
order accurate central-differences in the q and { directions. 

Results 

In order to validate the present nonequilibrium PNS 
code, two test cases were computed. The coordinate system 
employed in the present calculations is shown in Fig. 1. 

Test Case 1 

The first test case computed was that of hypersonic lam- 
inar flow of dissociating air over a 10 degree cone at 0 degree 
angle of attack. The altitude chosen was 60.96 km where 
the ambient pressure and temperature are 20.35 A’/m* and 
252.6 K , respectively. The remaining flow conditions are 

= 8100 m/s 

7*J, = 1200 K and noncatalytic wall 
d m - 0.2629 and c 3 „ = 0.7371 
Lt = 1.4 

The computation was started at x = 1.5 x 10~ s using 
the initial solution generated by the “stepback" procedure. 
The solution was then marched to x = 3.5. The march- 
ing Step site was chosen to be A£ = f6, where / > 1 and 
6, is the maximum thickness of the subsonic layer (in the 
present calculations, a value of 2 was assigned to /). The 
grid used in the calculations consisted of 67 points in the 
normal direction and 21 points in the meridional direction. 
The distance of the first point away from the body sur- 
face was varied linearly from 3 x 10" * to 2.3 x 10 -4 which 
determined the appropriate stretching parameter 0. The 
grid lines were placed normal to the body and the height 
of the outer boundary was kept fixed at 0.5. The edge of 
the boundary layer was located approximately using total 
enthalpy as the criterion. The edge values of pressure, tem- 
perature and velocity at the last station were then used as 
the uniform edge conditions for the reacting boundary-layer 
(RBL) code of Ref. 13. 
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(38) 


The surface pressure coefficient is defined as 


C„- 


Put P PC 
]p«>oc ! 


In Fig. 2, the axial variation of the surface pressure coeffi- 
cient is compared against the edge pressure of the boundary- 
layer code. The pressure predicted by the PNS code is ini- 
tially higher than the edge pressure of the RBI. code but 
eventually asymptotes to the latter value. The higher pres- 
sure is believed due to the leading edge effect. It must 
be recalled here that the PNS equations contain a normal 
momentum equation. This permits the interaction of the 
outer inviscid region with the inner viscous region. At the 
altitude considered, this interaction is fairly strong because 
the Reynolds number is low. Consequently, the pressure 
is higher. As the solution proceeds downstream, the in- 
teraction is reduced and one obtains near boundary-layer 
behavior. The increase in pressure leads to some differ- 
ences between the PNS and boundary-layer results as will 
be discussed later. 

The PNS equations are uniformly valid in the shock 
layer and as mentioned earlier the solution domain includes 
both the viscous and inviscid regions. In Figs. 3-8 only 20 
percent of the solution domain is shown in order to empha- 
size the details of the viscous boundary layer. Pitch plane 
profiles of tangential velocity and temperature at z ss 3.5 
are compared in Figs. 3 and 4. respectively. The agreement 
between the two codes is excellent. Mass fraction profiles 
of O and NO at x - 3.5 are compared in Figs. 5 and 6, 
respectively. The mass fractions have been normalised us- 
ing the corresponding wall values. The agreement between 
the two codes is again excellent. The wall values of the 
mass fractions were also found to be in very good agree- 
ment. The electron density (number of electrons per m*) 
is obtained from the mass fraction of the NO* ion. This 
follows from charge conservation and is defined as 


n: 


P' c no\ 
M NO ♦ 


(39) 


N 



Fig. 2 Wall pressure coefficient comparison. 




Temperature, T 


Fig. 4 Temperature profiles at x = 3.5. 


In Figs. 7 and 8, profiles of the mass fraction of NO~ 
and electron density obtained are compared. Again, the 
mass fractions and densities have been normalized using 
the corresponding wall values. The curve corresponding to 
the PNS calculations passes below the symbols representing 
the boundary-layer calculations. This is due to the fact that 
the PNS code predicts a higher amount of NO* near the 
wall. The onset of chemical reactions occurs earlier in the 
PNS code than in the RBL code because of the initial high 
pressures. This effect persists downstream and hence the 
disparity. The electron densities have been obtained from 
the mass fractions of NO* and consequently these densities 
will exhibit an identical behavior. 
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Fig. 8 Electron density profiles at x = 3.5. 



Fig. 6 SO mass fraction profiles at x = 3.5. 


The skin-friction coeflicienl and the Stanton number are 
defined as: 

<v v'i* W 

2 * <x> 


St - — 

~ HJ 


(«) 


where the wall shear stress ( Nfm ■*) is computed from 


r u = - P U 


av • ' 


■ Bn' I., 


(«) 


and the total heat transfer ( W/m J ) is computed from 


J . n . r— . . de, I 

** din ' I ' ^ 


*-\ 


(43) 


The first term in Eq. 43 is the conductive heating rate and 
the second the diffusive heating rate. The partial deriva- 
tive, djdn\ is taken in the direction normal to the surface 
of the body. In Figs. 6 and 10, the skin-friction coefficient 
and Stanton number are plotted as functions of the distance 
along the cone axis. The coefficients predicted by the two 
codes are in very good agreement. 
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Fig. 9 Skin-friction coefficient comparison. 



Fig. 10 Heat transfer coefficient comparison. 


The iota) drag coefficient is defined as 



(«) 


where D ’ (/V) is the total drag, i.e., the sum of the pres- 
sure drag and the skin-friction drag. The base pressure 
drag has been neglected. In the present study, trapesoidal 
integration was used to compute the total drag. The cross- 
sectional area at z = 3.5 was chosen as the reference area. 
The value of the drag coefficient predicted by the PNS code 
was 0.1037 which is in very good agreement (less than 1 per- 
cent) with the value of 0.1043 predicted by the RBL code. 


The ideal gas and equilibrium air models represent two 
physical extremes. In the former case, there are no reac- 
tions (infinitely slow reactions) and in the latter case, the 
reactions proceed at an infinite rate. In the gas model con- 
sidered in the present calculations, the reactions proceed at 
finite-rate and hence is between these two extremes. In or- 
der to demonstrate this, the PNS code of Ref. 20 was mod- 
ified and specialised to two-dimensional and axisymmetric 
bodies. New curve fits* 1 for the thermodynamic proper- 
ties of equilibrium air were used in the code. For the ideal 
gas calculations the ratio of specific heats was set to 1.4. 
The code was then used to compute the flow’ around the 
10 deg. cone for the same freestream conditions. In Fig. 
11, the temperature profiles at z = 3.5 obtained for the 
three gas models are compared. It is evident from the fig- 
ure that the peak temperature for the finite-rate chemistry 
case lies between the peaks for the ideal gas and equilib- 
rium air models. In the ideal gas case there are no internal 
degrees of freedom for the energy to be stored. For the 
finite-rate reaction case, some of the energy goes into excit- 
ing the internal degrees of freedom and some into chemical 
reactions. This effectively lowers the peak temperature. 
For the equilibrium air gas, thermochemical equilibrium is 
acheived instantaneously at every point of the flow field and 
hence a much lower peak temperature. In Fig. 12, the effect 
of the three gas models on the Stanton number is shown. 
The Stanton number has been plotted as a function of the 
axial distance. The equilibrium Stanton number is slightly 
greater than the finite-rate chemistry model which in turn 
is greater than the ideal gas Stanton number. 

Test Case II 

In order to validate the angle of attack capability of 
the code, a simple test case was chosen. In this test case, 
hypersonic laminar flow of reacting air over a 10 degree cone 



Fig. 11 Effect of gas models on temperature. 
Temperature profiles at z = 3.5. 
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Fig. 12 Effect of gas models on Stanton number. 
Axial variation of Stanton number. 


at 2.5 , 5, 7.5 and 10 degrees angle of attack was computed. 
The flow conditions for these calculations were the same as 
the flow conditions for Test Case 1. The starting solution for 
the PNS code for each angle of attack was obtained using 
the stepback procedure with the guessed initial solution 
being the frcestream. The mesh spacing on the windside 
had to be progressively refined for each angle of attack in 
order to properly resolve the boundary layer. The factor 
/ that multiplies the maximum subsonic layer thickness 
was decreased from 1.8 (2.5 deg case) to 1.25 (10 deg case) 
to maintain reasonable marching stepsizes. In each case 
the grid lines were placed orthogonal to the body. The 
starting solution was marched up to an axial location of z 
. = 2.5. Since there is a paucity of experimental data in such 
severe hypersonic regimes, the results of these calculations 
are simply compared against the results of the sero degree 
calculations. 

The computed results for this test case are too numer- 
ous to show in their entirety. Instead, only a representative 
sampling of the results is presented. In order to provide 
greater details of the flow field, only 50 percent of the solu- 
tion domain is shown in the figures that follow. The effect 
of angle of attack on the temperature is depicted in Figs. 
13a-.13d. The temperature profiles in the pitch plane of 
symmetry are shown in these figures. The following obser- 
vations can be made from these figures: (i) the boundary 
layer thickens considerably on the leeside while thinning on 
the windside. (ii) the edge temperalures increase rapidly 
on the windside while decreasing gradually on the leeside, 
and (iii) the peak temperatures decrease on the leeside but 
stay very nearly the same on the windside. The shock on 
the leeside weakens and begins to smear. At the highest 
angle of thp attack the shock cannot be discerned from the 
figure. This is due to a combination of central-differencing 
and coarseness of the grid. 


In Figs. Ha-14d, the effect of angle of attack on the 
mass fraction atomic oxygen in the pitch plane is shown. 
The amount of atomic oxygen at the wall decreases with 
increasing angle of attack. However, due to diffusional ef- 
fects atomic oxygen is present over a larger distance from 
the body. On the windside, the amount of atomic oxygen 
increases at the wall but is present over a smaller distance. 

For different angles of attack, the axial variation of the 
surface pressure and Stanton number on the windward and 
leeward meridians is shown in Figs. 15a-15d and Figs- 16a- 
16d, respectively. The increase in surface pressure and heat 
transfer on the windside and the decrease on the leeside 
with increasing angle of attack is evident from these figures. 

The computations were performed on either the CRAY- 
XMP/48 or the CRAY-2 (NAS) computers at NASA Ames 
Research Center. Each test case involved 10 x 10 block ma- 
trices and required 4.3 milliseconds per grid point per step 
on the CRAY-XMP computer. The first test case required 
2000 steps and the second test case required 1800 steps for 
the lowest angle of attack calculation to 600 steps for the 
highest angle of attack calculation. 

Concluding Remarks 

A new PNS code was developed to compute hypersonic 
laminar flow of chemically reacting air consisting of six 
species and electrons. The species were treated as calori- 
cally imperfect but thermally perfect gases. A noniterative, 
implicit, approximately-factored, space-marching method 
was used to solve the coupled set of gas dynamic and species 
conservation equations. In order to validate the code, two 
test cases were computed. The first test case was that of 
hypersonic flow over a 10 deg. cone at 0 deg. angle of at- 
tack. The results of this calculation were compared against 
those of a reacting boundary-layer code. The agreement 
was found to be very good. These results were further 
shown to lie between two extremes (ideal gas and equilib- 
rium air). The second test computed was that of hypersonic 
laminar flow of reacting air over a 10 deg. cone at several 
angles of attack. Due to the unavailability of experimental 
data for this case, the numerical results were simply com- 
pared against the results of the 0 deg. calculation. These 
calculations were performed not only to demonstrate the 
three-dimensional capabilities of the code but also to pro- 
vide “benchmark" calculations for future code developers. 
Work is currently being done in the area of shock fitting. 
A simple turbulence model (Cebeci-Smith) has also been 
incorporated into the present code and this option remains 
to be validated. The present study indicates that there is 
a need for a comprehensive experimental database for val- 
idation of nonequilibriuro codes. 
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(a) Fluxes 


Appendix A: Fluxes and Source Terms 


The (n + 4)-component inviscid flux vectors E 1 , P*, and G' are 


E' - {pu.pu* + p t puv,puu/,puff,puci,pucj,...,pue,,-i} T 
F* {pv, puv.pt; 2 + p,pvtv,pvff,pvci,pvc 2 ,. • .,pve B . i} r 
G' = {pw,puw,pvw,pu;*-t-p,pu;ff,pu;ei,ptt;ej,...,pu;e I ,-i} r 


(Al) 


14 


ORIGINAL PAGE IS 

OF, POOR QUALITY 


and the (n + d)-component viscous flux vector* B v , F*. and G* are 


B w = {0,r**,r*»,r“,ur** + vr“» + wr" - g*. mj.m},... ,m;_, } T 

F* = {0,T'’•,T•'',r»•,uT'’^ + vr" + vn»' - (* 2 ) 

G v = {0, r r**, r", ur** + *f ** + ter" - g*,mf,m5 mi-»} T 

The expression* for the component* of the (hear atresa tensor, the heat flux vector and the diffusion mass flux 
vector are given below 


2 

r** = -*i(2u, -v y - w„) 

ff* — - A*7i - - t»Jm; 

* 

* n = £/*(-«* + 2«>ir " w ») 

«» = -/*j*7V - £[h. - h n \m\ 
• 

2 

r" = - v y + 2u>.) 

g* = -UjxT, - - A.K 

r*» = *»(u r + v,) 

mj = fapD(c,), 

f“ = m(«» + »«) 

“ fitpDMu 

r** = p(», + u>„) 

"•I * 0»pP(e,) t 


The chemical source vector W' to 

W e = {0,0,0,0,0,tii,» J ,...,w, w«-i} r 


(A3) 


(Ad) 


where w, to the mass production or depletion rate of species *. 

(b) Chemical Sour ce Term* 

Consider a multicomponent system of n species undergoing m simultaneous elementary reactions. Let n t 
be the total number of reactants. These reactions can be represented symbolically as 

■» 

y > l,Ai s* 7 : v'kjAi *=l,2,...,m 

Is) im I 

where « the stoichiometric coefficients and A| is the chemical symbol of the Ith species. Using the 

law of mass action, the nondimensional mass production rate of species s to 


- •'UU/MT) ftw- - K>AT) ft|r»r|*- } (A5) 

km I ' rat rat ' 


The nondimensional mole-mass ratios of the reactants are defined as 


Tfr 


( e r /M, r= 1,2,. ...n 

E«sl ^(r-»)iiTf« f * B + 1,B + 2,...,tt| 


(d«) 


The reaction rates are functions of temperature and are expressed in the modified Arrhenius form, where 
the nondimensional forward and backward rates are written as 


K f , k [T) = exp(log,C;,* + -^+ C M log,D 
Kt,k(T) = expflog.O;,* + ^ + D,,*log.T) 


where 


CJ, ‘ = 

i>\,k - £(£?)* 


r r _ 

C7 }>4 _ __ 


T* 

1 oo 


n» - £!l* 

* oa 


(A7) 

(A8) 

(AO) 

(A10) 
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«* = II "l.r- A*£<r 


Mil) 


r*» I 


and Ci,*, Cj.*, Cj,*, /?*,*, and /?#,* are constants for a particular reaction i. 

(c) Transformed Fluxes 

Under the coordinate transformation, Eq. 18, the transformed fluxes and source vectors are given by 

ft* = (y)E' + (^)F' + (y)G' *' = + (^)F” + (^)G“ 

ft* = (y)R' + (^)F* + (y )G* ** = (y)E v + (^)G« 

6* = ( j)E' + (^)F* + ( j)G' = (j)B* + (J)F* + (^)G* 

W c = (i)W f 

J 


j _ 

~ 5(*,y,s) 


where J is the Jacobian of the coordinate transformation. 

The (n + 4)-component transformed viscous flux vector, f , for the thin-layer approximation is 

( 0 \ 

<1«» + Uv n + *»», 

<4«* + f»«e + 

<»«, + fe«e + <s««e 

i(ti - *»)(«*), + i(f, - 6)(u»), + $(*, - *»)(»’),+ 

d4(uv), + f»(utti),, + *e(t»»)e + tit f* + (*s - *?)£,!*« - An|(c»), 

*s( e i)e 
<e(**)e 


r-j 


where the coefficients £i,£*,...,<e are 


ti * P 
tt = P 
t»-n 


+ + (j)’ 

(yl’ + l’/l' + jl’y) 1 




**-s 

*-s 


(7><7> 


( jfXylJ 


M«) 


(A IS) 


(Alt) 


l * = l» = fopD 

,^ ) . + ( aL). + (t).] 


(d) Geometrical Parameters 

The Jacobian and the metrics of the coordinate transformation 

are 


J “ l*«(Vs*{ ~ Ms) + *n(Vt*t ~ VtU) + 1 

fc(V{*n - ‘ 

(A15) 

~J ~ (v**t " ~j - (*f«* _ x e*f) 

7 - (*sVf “ Me) 


7 = (Vf«t - v<*f) ^ = (*<»f - * f *<) 

7 = (*fV< - *m) 

(AW) 

7 = (Vt«e - Vn*<) 7 = ( x s*< - x <*e) 

7 = ( x <Ve - x eV<) 



Appendix B: Inviscid Jacoblans 

The inviscid Jacobian is an (n + 4) x (n + 4) matrix. The expressions for the elements of this matrix are 
given below. 

Define the following quantities 


fXC p ,7(l/.M. -!/*,,)-(*, 


- h n ) if s < n; 
if s = n. 


(Bl) 
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m 
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(*) Global Continuity 


(b) {-Momentum 


(c) gjjogwntum 



i. 

5, u + ^ y v + 5,w 


£» 


Ai,» *» 

Am = 5, 

a m = 5. 


A»,i * -uO + wx$«(y + A.) 
Aj,j = (1 - wx)5,u + 0 
Aaj = - wxi,® 

A*,« = £,u - u>x£,w 
Aj,» » u>x£* 


Aj, e = wx^i^i 

A*, r « U>x&tt a 


Aj,*+ 4 a uxi«^.-i 


(d) ^-Momentum 


(«) Energy 


A»,i = + wx^(y + **) 

Aj, a a S t v - uxS v u 
As^ = (1 - Wx)5 y V + # 

As.4 = £,V - U\£yW 

A*.ft a UxSy 


A*, ft a WX-S r dl 
Ag.r a WX^sdl 


A >( .44 » Wx5»d»-» 


* * . V* 

A 4 ,, = -wU + w X S,(y + d») 

A4,] a 5 r u; - wxS*«» 

A 4,3 = 5„VJ - lt)xS M V 


Aft.ft ~ WX^»^I 

A4,» = wx$,da 


A4,4 = (I - wx)5 a <l> 4 

Aft.ft a WXS* 



As,! a —0 H 
As,a = $mH 
Aft., » SyH 

Aft, 4 a 
Aft, ft = 0 


(f) Species Continuity (a a l,2,...,n - I) 


A(4«,i = -0e, 
* » 
Ai4«,a =yi*i 

Aft+«,» = ^|<| 

Aft+,,4 * ^|<» 

Aft+4,1+1 * ^ 
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Notes 

(1) The elements not listed above are identically sero. 

(2) The elements of the inviscid Jacobian dtf/dQ are obtained by setting 
St — (it S| = fy and S t — £,• 

(3) The elements of the inviscid Jacobian dt'/dQ are obtained by setting 
S t = »?«. S v = i)|, $| — 1 >», and w= 1. 

(4) The elements of the inviscid Jacobian dG'/dQ are obtained by setting 
St — (st = fyi ® it i end w - 1. 

Appendix C: Viscous Jacobian 

The viscous Jacobian is also an (n + 4) x (n + 4) matrix. The expressions for the elements of this matrix 
are given below. The following notation is used 

4>. = Jt. s = 1,2 8 (Cl) 

where It, it, ■..,(» are defined in Eq. AM. 

(a) (-Momentum 

Mj.i - ~ iV'i(^)s + ^(“)s + V>s(-)s) 

P P P 

; (C2) 

P 

Mr.< •- Vs(-), 

P 

(b) n-Momentum 


Ms,i = “ +lM~)e + ^•®(“)o] 

P P P 

Ms.s = 

Ms, a “ V»j(-)e 

P 

M 34 = ,M;)e 

P 

(c) (-Momentum 

Ms.i = + tM")e + ^s(“)s] 

M«, s - ^a(^)s 
P 

m 4 ,« = m'-U 

P 

(dj E^nergy 

Ms, i - ~\(4>t V’7)(--) n 4 (V'j - ^t){ — )* + (V>s - ^)("-)s + W«(— ),4- 

P P P P 

- 2<M— ), + M-U + (tfs - 0»)(A. - M(-)e+ 

p p p p 

(** " *t)(A» - *«)(-),+ ... + (*s - *r)(*»-. - *,)(—),] 

P P 

Ms.r = (V’i - V*r)(^), + V* 4 (-)„ 4 tM-)* 

P P P 

Ms,s = tM;)n 4 (V’J -• V>t)( V -)i» 4 tM-)e 

P P P 

Ms.s = V’st;), 4 0c(;), 4 (V»s - V»r)( 7 )s 
P P P 


(C3) 


(C4) 


(C5) 
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Ms,i — 0t(~ )i» 

P 

Ms, 8 * (th “ i ~ ^«)(")n 


Ms,»+4 = (<<>** 




fe) Specie* Continuity (a « 1,2,. ..,>»- 1) 


*•+..! « -*(T>e •-».* »“» 

P 

Ms+*,t+« = **= 1,2,. 1 

P 


(C6) 


Notes 

(1) The element* not lifted above ere identically sero. 

(*)()* = £(). 

Appendix D: Source Term Jacobian 

The chemical source vector, W*, is a function of the mixture temperature, density, and man concentrations. 
This is mathematically expressed as 


w« = jW‘(r,p7„p7» ry».) 


The Jacobian of the source term is 

. aw^ i aw* 
A = aq = J'aq 


Using Eq. Ol and the chain rule, the partial derivative in Eq. 02 can be written as 


aw* aw* ar aw* ar 
aqf “ ar aQ + ar ag 


(Dl) 

m 

m 


where T is the vector containing the man concentrations of the reactants. The derivative dW e /dT is easily 
evaluated since the reaction rate constants are the only quantities that depend explicitly on the temperature. 
The derivative dT/d Q is a I x (n 4 4) row vector whose elements are 


ar j f v* \ 

aQ = pc — \ 2 *’<■ u>,l,-(Aj - h H ),-(kt - h«),...,-(A B _i - A«)| (04) 

Evaluation of the second term of Eq. 03 is little more involved. The derivative dW*/dr is obtained by 
differentiating Eq. AS with respect to the elements of T and the derivative dT/dQ is easily evaluated using 
Eq. A6. 
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